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ABSTRACT 

We  present  a  methodology  for  the  efficient  calculation  of 
the  shock  Hugoniot  using  standard  molecular  simulation 
techniques.  The  method  is  an  extension  of  an  equation 
of  state  methodology  proposed  by  Erpenbeck  [J.  J. 
Erpenbeck,  Phys.  Rev.  A  46,  6406  (1992)]  and  is 
considered  as  an  alternative  to  other  methods  that 
generate  Hugoniot  properties.  We  illustrate  the 
methodology  for  shocked  liquid  N2  using  two  different 
simulation  methods:  (a)  the  Reaction  Ensemble  Monte 
Carlo  method  for  a  reactive  system;  and  (b)  the 
molecular  dynamics  method  for  a  non-reactive  system. 
The  method  is  shown  to  be  accurate,  stable  and  generally 
independent  of  the  algorithm  parameters.  We  find 
excellent  agreement  with  results  calculated  by  other 
previous  simulation  studies.  The  results  show  that  the 
methodology  provides  a  simulation  tool  capable  of 
determining  points  on  the  shock  Hugoniot  from  a  single 
simulation  in  an  efficient,  straightforward  manner. 


1.  INTRODUCTION 

The  behavior  of  materials  under  conditions  of  extreme 
temperature  and  pressure  is  of  significant  interest  in 
many  fields  of  physics  and  fluid  science  [1-4].  Of 
special  interest  are  energetic  materials,  a  class  of 
materials  of  critical  industrial  and  military  importance. 
These  materials  exhibit  chemically  and  physically 
interesting  behavior  when  exposed  to  extreme 
temperatures  and  pressures.  In  particular,  when 
subjected  to  shock,  energetic  materials  often  undergo 
rapid  reactions  that  produce  a  heterogeneous  mixture  of 
chemical  species  that  are  accompanied  by  huge  energy 
releases  and  can  produce  pressures  up  to  several  hundred 
GPa  and  temperatures  exceeding  10,000  K.  For  a 
sufficiently  strong  shock,  a  supersonic,  self-propagating 
reaction  wave  known  as  a  detonation  can  be  initiated. 
Unfortunately,  the  extreme  conditions  along  with  the 
short  time  and  length  scales  over  which  a  detonation 
occurs  poses  considerable  experimental  challenges  in 
characterizing  the  material  behind  the  detonation  front. 
Therefore,  a  concerted  effort  that  combines 
experimental,  theoretical,  and  simulation  approaches  is 
essential  for  furthering  our  understanding  of  such 
shocked  systems.  Advances  in  experimental  capabilities 
provide  us  with  crucial  property  data,  while  the 
continuing  development  of  accurate  equations  of  state 
have  allowed  reasonable  predictions  of  various  shock 


properties  [5,  6].  Similarly,  the  development  of  novel 
methods  to  simulate  these  complex  systems  has  been  the 
focus  of  research  efforts  and  has  recently  led  to  the 
invention  of  some  uniquely  effective  simulation  tools  [7- 
11].  These  classical  simulation  methods  can  be 
implemented  irrespective  of  rate  limitations,  the 
production  of  huge  energy  releases,  or  extreme 
thermodynamic  conditions. 

The  Hugoniot  curve,  a  commonly  calculated  property  in 
shock  and  detonation  science,  reveals  many  properties  of 
shocked  materials  and  knowledge  of  which  is  critical  to  the 
design  of  new  materials  and  application  platforms.  This 
curve  consists  of  the  set  of  ( PVT)  points  for  which  the 
Hugoniot  expression: 

Hg  =  E  -  E0  -  y2(P  +  P0)(F0  -V)  (1) 

is  zero.  In  Equation  (1),  E  is  the  specific  internal  energy,  P 
is  the  pressure,  and  V=l/p  is  the  specific  volume  (p  is  the 
specific  density).  The  term  specific  refers  to  the  quantity 
per  unit  mass,  while  the  subscript  “o”  refers  to  the  quantity 
in  the  initial  unshocked  state. 

Presently,  three  approaches  exist  for  calculating  the  shock 
Hugoniot  states  from  classical  molecular  simulation,  each 
with  their  own  advantages  and  disadvantages.  The  first 
approach,  which  we  term  here  the  Erpenbeck  equation  of 
state  method  (E-EOS),  is  the  most  indirect  of  the 
approaches.  The  original  version  of  the  method  involves 
performing  several  separate  simulations  at  appropriately 
chosen  temperatures  and  pressures.  Each  simulation 
generates  an  equation  of  state  point  for  subsequent 
evaluation  of  the  Hugoniot  expression  followed  by 
interpolation  to  locate  the  point  at  which  the  expression  is 
zero.  The  molecular  dynamics  (MD)  method  has  been 
implemented  in  the  Erpenbeck  approach  using  reactive 
potentials  that  mimic  chemical  bond  breaking  and  forming 
between  species  [12,  13]. 

An  approach  for  calculating  the  shock  Hugoniot  properties 
from  classical  molecular  simulation  that  is  more  direct 
than  the  E-EOS  method  is  the  piston-driven  molecular 
dynamics  method  [7,  8].  Piston-driven  MD  generates  a 
point  on  the  Hugoniot  curve  from  a  single  simulation,  thus 
avoiding  the  need  to  calculate  several  EOS  points  (as  in 
the  E-EOS  approach)  in  order  to  obtain  the  desired  result. 
The  method  mimics  the  laboratory  system  by  calculating 
properties  behind  the  shock  discontinuity  in  a  shock  wave 
simulation.  Shock  waves  are  produced  by  hitting  the  free 
edge  of  the  molecular  solid  with  a  rigid  layer  of  atoms  that 
are  moving  at  a  constant  velocity.  Different  shocked  states 
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are  obtained  by  starting  with  different  initial  piston 
velocities  [8]. 

A  third  approach,  termed  the  uniaxial  Hugoniostat 
method  [9],  is  a  molecular  dynamics  method  that  utilizes 
modified  equations  of  motion  that  constrain  the  system 
to  states  that  correspond  to  points  on  the  shock  Hugoniot 
curve.  This  method  is  computationally  more  efficient 
than  both  MD  implementations  of  the  E-EOS  method 
and  the  piston-driven  shock  wave  simulations.  The  E- 
EOS  method  requires  a  system  of  only  a  few  hundred 
atoms  (with  periodic  boundary  conditions  imposed),  but 
several  simulations  are  required  to  generate  a  single 
shock  Hugoniot  point.  The  piston-driven  MD  method 
will  produce  results  from  a  single  simulation,  however, 
the  system  size  must  be  sufficiently  large  so  that  the 
properties  behind  the  propagating  shock  wave  can  be 
averaged  over  time.  The  uniaxial  Hugoniostat  method, 
on  the  other  hand,  can  generate  a  point  on  the  Hugoniot 
curve  from  a  single  simulation  whose  system  size  is 
relatively  small. 

Unfortunately,  none  of  these  MD  methods  can  be  applied 
to  the  calculation  of  the  shock  Hugoniot  locus  over  a 
wide  range  of  temperatures  and  pressures  unless  certain 
conditions  are  fulfilled.  Most  energetic  materials 
respond  to  shock  by  decomposing  into  a  complex 
(sometimes  heterogeneous)  mixture  of  many  different 
chemical  species.  Thus,  for  these  multi-component 
systems,  the  implementation  of  the  uniaxial  Hugoniostat 
method,  the  MD  E-EOS  method  or  the  piston-driven  MD 
method  requires  either:  (1)  a  priori  knowledge  of  the 
relative  concentrations  of  each  chemical  species  in  the 
shocked  state;  or  (2)  a  reactive  potential  that  simulates 
bond  breaking  and  bond  formation.  Typically,  the 
relative  species  concentrations  of  the  shocked  state  are 
lacking,  moreover,  knowledge  of  these  quantities  is 
desired.  Furthermore,  although  significant  advances 
have  been  made  in  developing  reactive  potentials  for 
shocked  materials,  the  potentials  are  presently  limited  to 
idealized  representations  of  the  chemistry  that  occurs 
[< e.g .,  13].  The  most  accurate  interaction  potentials  for 
energetic  materials  available  at  this  time  are  all  non¬ 
reactive. 

Monte  Carlo  methodologies  circumvent  some  of  the 
restrictions  associated  with  the  MD  methods  for 
calculating  shock  Hugoniot  states.  The  Reaction 
Ensemble  Monte  Carlo  method  (RxMC)  [10]  and  the 
Composite  Monte  Carlo  method  [11]  have  both  been 
used  to  calculate  Hugoniot  properties  through  the  E-EOS 
approach.  These  two  closely-related  methods  do  not 
require  a  reactive  potential  or  a  priori  knowledge  of 
species  concentrations  for  each  Hugoniot  state.  They 
also  do  not  require  the  specification  of  species  chemical 
potentials  or  chemical  potential  differences  to  determine 
chemical  equilibrium  states  of  the  reactive  mixtures. 
Both  methods  have  been  applied  to  simple,  spherically- 
averaged  intermolecular  potentials  [10,  11]  but  can 


readily  be  applied  to  complex  potentials  that  include  multi¬ 
site  and/or  electrically-charged  species  as  well  as  multi¬ 
phase  mixtures.  Therefore,  in  the  absence  of  reactive 
potentials  or  a  priori  knowledge  of  species  concentrations, 
the  only  applicable  approaches  for  simulating  the  Hugoniot 
properties  of  a  shocked  material  are  the  Erpenbeck  EOS 
method  performed  using  either  the  RxMC  or  Composite 
MC  methods. 

As  previously  mentioned  however,  the  original  E-EOS 
method  requires  simulations  of  several  equation  of  state 
points  to  generate  a  single  point  on  the  shock  Hugoniot 
curve.  Each  separate  simulation  requires  sufficient 
equilibration  and  data  collection  steps.  In  an  effort  to 
reduce  the  number  of  steps  and  to  minimize  associated 
computational  costs,  we  have  implemented  a  numerical 
approach  within  the  framework  of  the  E-EOS  approach. 
The  resulting  method  requires  only  a  single  simulation  to 
determine  a  point  on  the  Hugoniot  curve.  The  fitting 
procedure  used  to  determine  the  root  of  the  Hugoniot 
expression  (Hg=0)  in  the  original  version  of  the  E-EOS 
approach  is  replaced  by  an  iterative  numerical  procedure 
built  into  the  framework  of  the  simulation. 

A  brief  illustration  of  the  method  using  isothermal-isobaric 
ensemble  ( NPT)  simulations  is  given.  The  simulation  is 
initiated  at  the  specified  temperature  and  pressure,  and  the 
Hugoniot  expression  is  evaluated  (using  instantaneous 
values  of  V,  P  and  E  that  depend  only  on  the  current 
configuration)  at  periodic  intervals  during  the  simulation 
run  and  accumulated  for  averaging.  The  averaged 
Hugoniot  value  is  then  used  in  a  numerical  root-finding 
algorithm  (e.g.,  Newton-Raphson  [14])  to  provide  an 
estimate  of  the  Hugoniot  pressure.  The  imposed  pressure 
for  the  simulation  is  subsequently  changed  to  correspond 
to  the  new  estimate  of  the  Hugoniot  pressure.  The 
simulation  continues  using  this  new  imposed  pressure 
constraint.  This  process  is  repeated  until  the  Hugoniot 
function  converges  to  zero  (more  precisely,  within  a 
desired  tolerance  of  zero).  The  value  of  the  pressure  and 
corresponding  volume  averaged  over  the  entire  simulation 
run  characterize  the  Hugoniot  state  at  that  temperature. 

The  method  is  akin  to  the  phase  equilibria  methods  that 
utilize  thermodynamic  integration  to  determine 
coexistence  behavior  [15-17].  In  these  calculations,  a 
finite-difference  algorithm  is  used  to  numerically  integrate 
the  differential  equations,  the  Clausius-Clapeyron  [15,  16] 
or  the  Gibbs-Duhem  expressions  [17],  which  govern  the 
changes  in  thermodynamic  parameters  along  the  phase 
coexistence  curve.  Similarly,  for  the  method  introduced 
here,  a  numerical  estimate  of  the  root  of  the  Hugoniot 
expression  is  made.  Both  approaches  are  iterated  until  the 
desired  convergence  has  been  reached. 

In  summary,  we  demonstrate  the  accuracy  and  stability  of 
a  method  to  calculate  the  shock  Hugoniot  properties  of 
materials  when  implementing  the  E-EOS  method  either  in 
an  MD  or  RxMC  framework.  (Implementation  of  the 


modified  E-EOS  method  using  the  Composite  MC 
method  is  analogous  to  the  RxMC  method  and  will  not 
be  demonstrated  here.)  The  method,  which  we  term  the 
adaptive  Erpenbeck  equation  of  state  method  (AE-EOS), 
is  intended  to  be  a  tool  that  is  an  alternative  to  the 
existing  methodologies  in  order  to  overcome  some  of 
their  limitations.  We  demonstrate  the  validity  of  the  AE- 
EOS  method  for  calculating  the  Hugoniot  properties  of 
liquid  N2  in  the  non-reactive  regime  using  molecular 
dynamics  and  in  the  reactive  regime  using  Reaction 
Ensemble  Monte  Carlo.  The  outline  of  the  paper  is  as 
follows.  The  formalism  and  practical  details  of  the 
methodology  are  presented  in  Section  II.  Applications  of 
the  method  are  given  in  Section  III,  while  assessments  of 
the  results  are  given  in  Section  IV. 


II.  METHODOLOGY 


A.  Formalism 


The  Hugoniot  function  Hg=Hg[2s,  7,  P  (or  V)],  and  so  a 
search  for  Hg=0  in  the  original  Erpenbeck  EOS  method  is 
implemented  by  fixing  two  independent  variables  in  a 
simulation  and  calculating  the  remaining  variable.  For 
example,  when  simulating  with  an  isothermal-isobaric 
ensemble,  E  is  calculated  from  the  simulation. 
Typically,  several  simulations  must  be  performed  using 
various  choices  of  the  independent  variables  to  generate 
sufficient  points  so  that  the  Hugoniot  state  can  be 
obtained  through  interpolation.  (Detailed  outlines  of  the 
original  E-EOS  approach  using  the  molecular  dynamics 
and  the  RxMC  methods  can  be  found  in  Refs.  [12]  and 
[10],  respectively.)  The  adaptive  Erpenbeck  equation  of 
state  method  presented  in  this  work  eliminates  the 
interpolation  procedure  by  using  a  root- finding  algorithm 
to  determine  Hg=0.  For  example,  the  expression  for 
finding  the  root  of  a  function  using  the  Newton-Raphson 
method  is  [14] 


n+ 1 


/(*J 

f'{Xn) 


(2) 


where/' (xj  =  df(xn)/dxn  ,Axn)= Hg  and  we  choose  either 

xn=E ,  T,  or  P(or  V).  The  basis  of  the  AE-EOS  method  is 
to  begin  a  simulation  at  an  initial  xm  and  after  a 
prescribed  number  of  simulation  steps,  the  quantities 
J[xn)  and f\xn)  are  calculated  and  xn+1  is  determined.  The 
simulation  continues  using  the  predicted  value  xn+1.  This 
procedure  is  repeated  until  the  results  converge. 

In  the  following,  we  demonstrate  the  AE-EOS  method 
for  P  as  the  independent  variable  used  to  find  the  root  of 
the  Hugoniot  function.  The  working  expression  for 
predicting  the  pressure  of  the  Hugoniot  state  within  the 
framework  of  the  Newton-Raphson  procedure  [Eq.  (2)] 
is  then 


ppredicted  _  pcurrent 


(3) 


In  the  procedure  presented  here,  the  predicted  pressure, 
ppredicted  ?  [ s  determined  using  averaged  instantaneous 
values  of  the  Hugoniot  expression  [Eq.  (1)]  and  its 
derivative  with  respect  to  pressure,  i.e.,  nc™rrent= 

and  dRcrnt_/dHg 
dP  \  dP 


where  ‘<  >’  denotes  averages  of  instantaneous  values.  The 
instantaneous  Hugoniot  values  (and  derivatives)  were 
determined  using  the  corresponding  instantaneous  values 
of  7,  E ,  and  T  generated  during  the  simulation. 

Next,  since  the  internal  energy  can  be  written  as  [10] 

Cf 

£  =  £  y ,H°  +  U£onf  -  RT  >  (4) 

1=1 


then  Eq.  (1)  can  be  rewritten  as 


H*  = 


(( c>  \  i 

Xy,7/,0+t/COTf-Rr  -E0  -±{PV0-PV  +  P0V0-VP0 

V  V /=1  )  ) 

(5) 


where  E0,  P0,  and  V0  are  the  values  of  the  initial  state  (and 
thus  constant)  and  ^conf  'YjJv{ri.)  w^ere  Uy  is  the  pair 

i  j>i 

potential  energy  [18].  All  quantities  in  Eqs.  (4)  and  (5)  are 
used  here  in  the  context  of  instantaneous  quantities  that 
depend  only  on  the  current  configuration.  The  derivative 
term  required  in  Eq.  (3)  is  then 


H - — —  (u  conf  )  — 

dP  ’ 


d_ 

dP 


(PVo-PV  +  P0V0-VP0 


(6) 


The  first  and  third  terms  on  the  r.h.s.  of  Eq.  (6)  are  not 
functions  of  P  and  can  be  eliminated.  Further  since  lfonf  is 
calculated  as  an  instantaneous  value,  <[7conf=iyconf(r)  only 
and  thus  can  be  eliminated.  Finally,  the  last  term  is  readily 
solved,  so  that  Eq.  (6)  reduces  to 


dH  ,  , 

- !L  =  UV-V) 

dP  2V  °; 


(7) 


The  algorithm  for  the  adaptive  Erpenbeck  equation  of  state 
method  with  P  chosen  as  the  independent  variable  is  as 
follows: 


Step  1:  Set  the  temperature  for  the  Hugoniot  state  (7Hg). 
Step  2:  Guess  the  pressure  for  this  Hugoniot  state 

^ pcurrent ^ 

Step  3:  Perform  an  isothermal-isobaric  ensemble 
simulation  (MD  or  RxMC)  at  7Hg  and  pcurrentm 


Step  4:  After  allowing  the  system  to  relax  to  pcurren\ 
accumulate  instantaneous  values  of  Hg  and 
dRg/dP  during  the  simulation  using  Eqs.  (5)  and 
(7),  respectively. 

Step  5:  After  a  prescribed  number  of  steps,  calculate 
averaged  values  of  Hg  and  dRg/dP  and  predict 
the  Hugoniot  pressure  using  Eq.  (3).  (The 
averaged  values  of  Hg  and  dRg/dP  can  be 
determined  by  several  methods,  which  are 
considered  in  the  next  section.) 

Step  6:  Repeat  steps  (3)-(6)  until  the  results  converge  to 
the  desired  statistical  uncertainty. 


B.  Practical  Details 

Next,  we  consider  a  few  practical  details  of 
implementing  the  AE-EOS  method.  Our  intent  is  to 
generalize  the  method  for  implementation  into  any  of  the 
standard  molecular  simulation  techniques  (MD,  RxMC, 
or  Composite  MC).  We  consider  the  effect  of  several 
parameters  on  the  accuracy  and  stability  of  the  method. 
Below  we  address  these  issues,  the  logic  behind  our 
choices,  and  the  tradeoffs  involved. 

There  exists  a  wide  range  of  root-finding  algorithms, 
including  the  Newton-Raphson  method,  the  Secant 
method,  the  Bi-section  method,  and  Halley’s  method 
[14].  Newton-Raphson  (Eq.  (2))  is  a  rather 
straightforward  method  but  can  be  unstable  near  a 
horizontal  asymptote  or  local  minimum.  A  similar 
algorithm  is  Halley’s  method  that  includes  an  additional 
term  ( df'(xn)/dxn )  from  the  Taylor  series  in  the 

derivation  of  the  method.  When  the  pressure  is  chosen 
as  the  independent  variable,  <72Hg/<7P2=0  (see  Eq.  (7))  so 
Halley’s  method  reduces  to  the  Newton-Raphson 
method.  Another  root  finding  method  is  the  Secant 
method,  which  estimates  the  derivative  term  using 
^  _  f(xn)  ~f(xn-\) .  However,  use  of  the  method  in 
xn  -  xn_} 

the  AE-EOS  method  requires  two  recent  points  along  the 
Hugoniot  curve  as  opposed  to  only  one  for  the  Newton- 
Raphson  method.  Finally,  if  we  can  be  certain  that  the 
solution  of  the  Hugoniot  expression  lies  within  a  known 
interval,  then  we  can  iteratively  converge  to  the  solution 
using  the  Bi-section  method.  However,  a  balance  must 
be  established  between  statistical  uncertainty  and  the 
desired  convergence  when  implementing  the  Bi-Section 
method  in  a  molecular  simulation,  since  statistical 
fluctuations  cannot  be  greater  than  the  size  of  the 
interval.  In  an  effort  to  keep  the  method  proposed  here 
as  general  and  straightforward  as  possible,  we  have 
implemented  the  Newton-Raphson  method.  The  well- 
known  problems  of  this  method  near  a  local  minimum  or 
asymptote  have  not  been  encountered  for  the  Hugoniot 
expression  in  this  work  as  well  as  for  other  work  [10,  12, 
13],  but  one  should  be  mindful  of  its  limitations. 


In  any  molecular  simulation,  it  is  necessary  to  design  a 
starting  configuration  so  that  the  relaxed  system  is 
physically  reasonable  and  computationally  consistent.  For 
example,  consider  an  isothermal-isobaric  ensemble 
simulation  where  periodic  boundaries  are  imposed  and 
where  the  potential  energy  function  has  a  limited 
interaction  range.  In  such  a  case,  an  appropriate  number  of 
molecules  must  be  chosen  so  that  the  relaxed  box  size  is 
consistent  with  the  minimum  image  convention,  i.e.,  one- 
half  the  box  size  must  be  greater  than  or  equal  to  the 
potential  cutoff  distance  [18].  Similarly,  appropriate 
starting  conditions  for  the  AE-EOS  method  are  required, 
particularly  for  the  initial  guess  of  the  imposed  pressure. 
Since  the  converged  result  (i.e.  Hg=0)  will  produce  a 
Hugoniot  pressure  that  is  equal  to  the  imposed  pressure 
(within  some  specified  tolerance),  it  is  desirable  to  choose 
an  initial  pressure  that  is  a  good  estimate  of  the  actual 
Hugoniot  pressure.  Although  we  will  show  in  Section  III 
that  the  AE-EOS  method  is  largely  insensitive  to  the  initial 
pressure  guess,  extremely  poor  initial  guesses  could  result 
in  numerical  failures.  For  example,  consider  the 
simulation  of  a  Hugoniot  state  such  that  THg>T0  (recall  that 
rHg  is  the  temperature  of  the  chosen  Hugoniot  state  and  T0 
is  the  temperature  of  the  unshocked  material).  If  the  initial 
guess  of  the  pressure,  Pinitiai  guess?  causes  an  expansion  of  the 
simulation  cell  such  that  the  specific  volume  is  larger  than 
the  specific  volume  of  the  unshocked  material,  a  negative 
pressure  value  will  be  predicted.  From  a  practical 
standpoint,  this  is  an  unphysical  occurrence  since  it  implies 
that  the  material  has  expanded  upon  shock  rather  than 
being  compressed.  Furthermore,  from  a  computational 
standpoint,  the  simulation  cell  will  never  converge  to  a 
negative  imposed  pressure.  Such  an  occurrence,  however, 
is  analogous  to  choosing  a  starting  configuration  that 
relaxes  to  a  physically  unreasonable  and  computationally 
inconsistent  state. 

Consider  the  following  ad  hoc  approach  to  choosing  a 
reasonable  initial  guess  of  the  pressure  for  Eq.  (3).  First, 
assume  that  the  shocked  material  does  not  decompose  (i.e., 
chemically  react).  This  is  a  reasonable  approximation  at 
low  shock  pressures  and  reduces  the  first  term  on  the  r.h.s 
of  Eq.  (5)  to  PP starting  material-  Next,  neglect  the  contribution 
of  I/onf  and  thus  eliminate  the  second  term  on  the  r.h.s  of 
Eq.  (5).  This  approximation  has  no  physical  justification, 
however,  a  short  simulation  could  be  performed  to 
calculate  lfonf  (although  probably  unnecessary  given  the 
lack  of  sensitivity  of  the  final  result  on  the  initial  guess  of 
the  pressure).  Finally,  estimate  the  amount  of  compression 
the  starting  material  will  undergo,  e.g.,  V=0.7Vo.  This 
estimate  presumably  can  be  predicated  on  previous  studies 
of  the  material  or  similar  materials.  With  these 
approximations,  an  initial  estimate  of  the  Hugoniot 
pressure  can  be  determined.  Furthermore,  as  points  along 
the  Hugoniot  curve  are  determined,  better  estimates  of  the 
initial  pressure  can  be  made  by  using  these  Hugoniot 
states. 


For  completeness  of  study,  we  consider  three  different 
schemes  for  averaging  the  instantaneous  values  of  Hg  and 
dRg/dP ,  which  are  then  used  in  Eq.  (3):  (1)  block 
averages;  (2)  running  averages;  and  (3)  block-to-running 
averages.  Block  averages  are  taken  from  a  limited 
number  of  configurations  immediately  preceding  the 
pressure  adjustment  step,  while  running  averages  are 
taken  continuously  over  all  configurations  generated 
during  the  simulation  run.  The  block-to-running 
averages  scheme  uses  a  block-averaging  scheme  for  the 
equilibration  period  of  the  simulation  run  and  then 
continues  with  a  running  average  scheme  for  the 
production  period.  This  scheme  may  most  effectively 
remove  the  effects  of  a  poor  initial  guess,  while  we 
expect  the  running  average  scheme  to  be  the  most 
effective  alternative  since  fluctuations  in  the  pressure 
will  be  become  increasingly  damped  as  the  simulation 
proceeds.  Block  averaging  methods  will  likely  be  more 
slowly  converging  at  best,  and  unstable  at  worst. 
Moreover,  running  average  schemes  have  been  the  most 
successful  scheme  in  the  finite-difference  algorithms 
used  in  the  phase  coexistence  methods  mentioned 
previously  [15-17]. 

We  also  consider  the  effect  of  the  frequency  of  re- setting 
the  pressure  during  the  simulation.  Less  frequent 
updates  are  expected  to  cause  the  results  to  converge 
more  slowly  while  more  frequent  updates  could  possibly 
cause  the  root-finding  scheme  to  become  unstable  or  to 
fluctuate  too  greatly. 

A  final  note  on  the  convergence  of  the  system  to  the 
predicted  pressure  is  considered.  Step  (4)  in  the 
algorithm  outline  allows  the  system  to  converge  to  the 
predicted  pressure  value  (within  a  few  %  of  the  predicted 
pressure  for  the  most  recent  simulation  steps)  before  re¬ 
evaluating  the  Hugoniot  expression  and  it’s  derivative 
(dRg/dP)  in  the  equilibration  period  only.  This  ensures 
that  even  for  large  changes  in  the  predicted  pressure, 
equilibrated  information  is  still  used  in  the  Hg  and 
dRg/dP  calculation.  Typically,  these  large  changes  will 
only  occur  during  the  earliest  stages  of  the  simulation. 
At  later  times  during  the  production  cycles,  this  criterion 
is  nearly  always  satisfied. 


III.  APPLICATION 

For  demonstrative  purposes,  several  shock  Hugoniot 
states  of  liquid  nitrogen  are  considered.  The  shock 
Hugoniot  properties  were  predicted  based  on  the  initial 
states  calculated  previously:  £=77.0  K;  p=0.808  g/cm3; 
£=50.49  MPa;  £=-0.441  kJ/g  [10].  At  pressures  higher 
than  ~30  GPa  along  the  Hugoniot  curve,  the  dissociation 
reaction  of  molecular  nitrogen  (N202N)  occurs. 
Therefore,  we  demonstrate  the  AE-EOS  method  using 
the  molecular  dynamics  technique  only  at  pressures 
below  30  GPa  while  we  demonstrate  the  AE-EOS 
method  using  the  RxMC  method  at  a  wider  range  of 


pressures.  Particles  interact  through  an  exponential-six 
potential,  where  potential  parameters  are  given  in  Ref. 
[19].  A  spherical  cutoff  for  the  particle-particle 
interactions  was  applied  at  2.5 rW;N2  with  long-range 
corrections  added  to  account  for  interactions  beyond  this 
distance  [20].  Electrostatic  interactions  between  species 
were  ignored.  The  unlike  interactions  between  species  i 
and  j  were  approximated  by  the  Lorentz-Berthelot 
combining  rules  [21].  3375  N2  molecules  were  used  with 
all  calculated  quantities  reduced  by  the  exponential-6 
potential  energy  (d)  and  size  (rm)  parameters  of  N2. 
Periodic  boundary  conditions  were  imposed  for  all 
dimensions.  Thermochemical  reference  data  were  used  in 
calculating  the  ideal-gas  enthalpies  (//z°)  required  in  Eq. 
(5)  [22,  23]. 


A.  Molecular  dynamics 

Molecular  dynamics  simulations  in  the  isothermal-isobaric 
ensemble  were  performed  using  the  Leap-Frog  Verlet 
algorithm  [18,  20]  and  the  Melchionna  modification  of  the 
Hoover-Nose  equations  of  motion  [24].  A  thermostatting 
rate  of  50  ps"1  was  used  to  maintain  the  imposed 
temperature  while  a  barostatting  rate  ranging  from  0.032- 
0.042  ps"1  was  used  to  maintain  the  imposed  pressure. 
Initial  configurations  were  generated  from  a  face-centered- 
cubic  (fee)  lattice  structure  with  initial  particle  velocities 
selected  from  a  Boltzmann  distribution  that  corresponded 
to  the  imposed  temperature.  Preceded  by  an  equilibration 
period  of  0.127-0.254  ns  during  which  the  pressure  was 
not  re-set  using  Eq.  (3),  trajectories  were  followed  for  1.32 
ns  with  time  steps  ranging  from  0.00763-0.0102  ps.  All 
pressure  values  reported  were  determined  using  the  virial 
theorem  [18]. 

Three  state  points  along  the  shock  Hugoniot  curve  were 
determined:  £=883.9;  3912.4;  6778.1  K.  These  state 
points  are  below  the  regime  in  which  N2  dissociates  into 
atomic  nitrogen.  For  each  point,  the  Hugoniot  pressure 
(£ Hg)  was  predicted  in  two  simulations,  one  of  which  the 
initial  pressure  was  much  lower  than  the  Hugoniot 
pressure  and  one  in  which  the  initial  pressure  was  too  high. 
The  effect  of  the  frequency  of  re-setting  the  imposed 
pressure  was  also  studied.  Two  cases  were  considered,  re¬ 
setting  at:  (a)  every  100  steps;  and  (b)  every  500  steps. 
Following  the  initial  equilibration  period  used  to  relax  the 
system  from  the  fee  crystal  to  the  imposed  thermodynamic 
condition,  an  additional  0.305  ns  of  the  trajectory  was  used 
to  further  equilibrate  the  system  after  the  AE-EOS 
algorithm  is  implemented  ( i.e .  the  pressure  is  re-set  at 
specified  intervals).  All  quantities  calculated  during  this 
time  interval  were  not  included  in  the  final  averages.  A 
tolerance  value  of  +1-5%  was  used  in  Step  (4)  for  the 
pressure  (see  Section  II. A),  i.e.,  the  calculated  pressure 
was  required  to  be  within  +1-5%  of  the  most  recent  ^Hg 
prediction  before  re-evaluating  Hg  and  dRg/dP  and  re¬ 
setting  of  the  imposed  pressure. 


Table  I:  Predicted  shock  Hugoniot  states  of  liquid  N2  using  molecular  dynamics  in  the  AE-EOS  method.1 


B.  Reaction  Ensemble  Monte  Carlo 

The  Reaction  Ensemble  Monte  Carlo  method  was  used 
to  assess  the  accuracy  of  the  AE-EOS  method  at  a  wider 
range  of  conditions  than  considered  using  the  molecular 
dynamics  technique  including  conditions  under  which  N2 
dissociates  (N2<=>2N).  Details  of  the  methodology  can 
be  found  in  the  original  papers  [25-27]  as  well  as  in 
recent  applications  of  the  technique,  which  implemented 
the  E-EOS  method  [10].  In  addition  to  intermolecular 
potentials  that  describe  non-reactive  interactions  between 
species  N2  and  N  in  the  equilibrium  mixture,  RxMC  also 
requires  inputting  the  ideal-gas  internal  modes  (vibration, 
rotation,  electronic).  The  vibrational  and  rotational 
contributions  to  the  ideal-gas  partition  functions  of  both 
species  were  calculated  using  a  standard  source  [28],  and 
supplemented  with  electronic  level  constants  that 
included  the  ground  state  and  six  excited  electronic  states 
for  N2  [22],  along  with  electronic  energy  levels  for  N 
taken  from  Moore  and  Gallagher  [23]. 

Constant-pressure  RxMC  simulations  of  shocked  N2 
were  initiated  from  3375  N2  particles  placed  on  an  fee- 
lattice  structure.  Simulations  were  performed  in  steps, 
where  a  step  (chosen  with  equal  probability)  was  either  a 
particle  displacement,  a  forward  reaction  step,  or  a 
reverse  reaction  step.  A  change  in  the  simulation  cell 
volume  was  attempted  every  500  steps.  Simulations 
were  equilibrated  for  1.5xl06  steps  after  which  averages 
of  the  quantities  were  taken  over  8xl06  steps. 
Uncertainties  were  estimated  using  the  method  of  block 
averages  by  dividing  the  production  run  into  10  equal 
blocks  [20].  Reported  uncertainties  are  one  standard 
deviation  of  the  block  averages.  (Note  that  these  block 
averages  do  not  correspond  to  the  averaging  used  in  the 
^Hg  prediction  scheme.)  The  maximum  displacement 
and  volume  change  were  adjusted  to  achieve  an 
acceptance  fraction  of  approximately  0.33  and  0.5, 
respectively.  Depending  on  the  system  conditions,  the 
acceptance  fraction  of  the  reaction  steps  ranged  from 
0.075-0.375. 

Three  points  along  the  shock  Hugoniot  curve  were 
considered:  7=2008.4;  7963.03;  10935.2  K.  Again  for 
each  point,  .Png  was  predicted  from  two  initial  guesses 
that  differ  from  the  known  value  [10]  by  +/-  67%.  We 
predicted  the  shock  Hugoniot  properties  based  on  the 
same  initial  state  used  previously  [10]  and  in  the 
molecular  dynamics  study  above.  A  tolerance  value  of 
+/- 2.5%  was  used  in  Step  (4)  for  the  pressure  (see 
Section  II. A),  a  value  slightly  lower  than  used  in  the 
molecular  dynamics  simulations  above.  Analogous  to 
the  study  above,  the  effect  of  the  frequency  of  re¬ 
evaluating  the  Hugoniot  pressure  and  subsequent  re¬ 
setting  of  the  imposed  pressure  was  studied.  Two  cases 
were  considered,  re-setting  at:  (a)  every  5000  steps;  and 
(b)  every  50,000  steps.  Additionally,  three  different 
averaging  schemes  of  the  instantaneous  values  of  Hg  and 


dRg/dP  were  assessed:  (a)  block  averages;  (b)  running 
averages;  and  (c)  block-to-running  averages.  A 
description  of  each  type  is  given  in  Section  II.B.  For  each 
of  the  initial  pressure  guesses,  all  six  series  (2  Hg 
frequency  evaluations  X  3  averaging  schemes)  were 
considered. 

Comparisons  between  the  Hugoniot  properties  predicted 
by  the  E-EOS  and  the  AE-EOS  methods  for  the  7963.03  K 
case  are  given  in  Table  II.  We  find  excellent  agreement 
between  the  E-EOS  and  the  AE-EOS  results  for  all 
quantities  calculated.  Nearly  all  AE-EOS  quantities  fall 
within  +/-0.5%  of  the  E-EOS  calculations  with  none 
greater  than  +/-0.9%.  Consequently,  no  dependence  on  the 
initial  pressure  guess  or  the  Hugoniot  expression 
evaluation  frequency  is  evident.  Likewise,  no  apparent 
dependence  on  each  of  the  averaging  schemes  is  found. 


IV.  DISCUSSION 

We  have  presented  a  computationally  efficient 
methodology  for  calculating  the  shock  Hugoniot  properties 
of  materials  using  classical  molecular  simulations.  The 
method  is  an  extension  of  the  Erpenbeck  EOS  approach 
and  allows  for  the  determination  of  a  point  along  the 
Hugoniot  curve  from  a  single  simulation.  The  method, 
termed  the  adaptive  Erpenbeck  equation  of  state  method 
(AE-EOS),  uses  a  numerical  estimate  of  the  root  of  the 
Hugoniot  expression  to  determine  the  corresponding 
thermodynamic  state.  AE-EOS  is  applicable  for  any 
simulation  method  that  determines  points  along  the 
Hugoniot  curve  by  generating  EOS  points,  which  includes 
molecular  dynamics  [18,20],  Reaction  Ensemble  Monte 
Carlo  [25,26],  and  the  Composite  Monte  Carlo  method 
[11].  The  method  was  demonstrated  to  be  accurate  and 
numerically  stable.  For  the  systems  in  this  study,  the 
method  was  not  particularly  sensitive  to  the  algorithm 
parameters  indicating  the  method’s  robustness  and  ability 
to  be  readily  implemented.  Furthermore,  a  substantial 
savings  in  the  computational  requirements  is  gained 
through  the  implementation  of  the  AE-EOS  method.  For 
comparable  statistical  uncertainty,  4-5  complete 
simulations  are  required  to  generate  a  single  point  on  the 
Hugoniot  curve  for  the  original  Erpenbeck  method  while 
only  a  single  simulation  is  needed  in  the  AE-EOS  method. 
Hence,  a  computational  gain  of  approximately  4-5  fold  is 
made  when  implementing  the  AE-EOS  method.  Overall, 
the  AE-EOS  method  appears  to  be  slightly  more  accurate 
and  stable  when  applied  to  the  RxMC  method  as  compared 
to  the  molecular  dynamics  method.  This  may  be  attributed 
to  the  relative  stability  of  the  Monte  Carlo  method  over  the 
molecular  dynamics  technique  for  constant-temperature 
and  constant-pressure  simulations. 

The  AE-EOS  method  is  intended  as  an  alternative  tool  to 
similar  techniques,  namely,  piston-driven  molecular 
dynamics  [7,8]  and  uniaxial  Hugoniostat  molecular 


dynamics  [9],  since  AE-EOS  offers  some  notable 
advantages  when  simulating  reactive  mixtures.  The  most 
striking  advantage  is  that  the  AE-EOS  method  does  not 
require  either:  (1)  a  priori  knowledge  of  the  relative 
concentrations  of  each  chemical  species  in  the  shocked 
state;  or  (2)  a  reactive  potential  that  simulates  bond 
breaking  and  bond  formation.  At  least  one  of  these 
conditions  must  be  met  for  the  other  methods  to  be 
applicable  to  the  simulation  of  reactive  mixtures. 
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Motivation 


Need  for: 

1 .  Better  performing  energetic  materials 

2.  Materials  to  keep  pace  with  advancing  technology 

(notional  materials) 


Goal: 


Efficient  design  of  energetic  materials 


Design  Approaches 

1 .  Laboratory  experiments 

•  Trial  and  error  process 

•  Years  to  develop  ~  20  years  for  M43  propellant 

•  Costly  &  environmentally  50,000  lbs  XM39  propellant 

harmful  (Millet  et  al  Clean  Prods.  &  Proc.  2000) 

2.  Theoretical  models  (Ex:  Cheetah) 

•  Fast  calculations 

•  Limited  predictive  capabilities  require  EOS 

•  Not  practical  for  notional  materials 

3.  Computational  methods 


Computational  Methodologies 


Based  on  SDSC  Blue  Horizon  (SP3) 
612-1024  processors 
1.728  Tnops  peak  performance 
CPU  time  »  1  week  /  processor 


Continuum 

Methods 
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CFD 

Cyclops 


Molecular  Simulation 


Predict  macroscopic 
properties  from 
molecular-level  behavior 


►  Simulation  of  chemically  reacting  systems: 


1 .  Molecular  Dynamics : 

v'  Requires  reactive  potential  for  each  reaction 
S  Not  applicable  to  mixtures 


2.  Reaction  Ensemble  Monte  Carlo 


Reaction  Ensemble  Monte  Carlo 


►  simulates  chemical  reaction  equilibria 


►  not  limited  by  reaction  rate 


►  not  limited  by  activation  energy  barrier 


❖  J.K.  Johnson,  K.  Gubbins,  and  A.  Panagiotopoulos  (1994) 

❖  W.  Smith  and  B.  Triska  (1994) 


Reaction  Ensemble  Monte  Carlo 


ideal  gas  reaction 


RxMC 


predict  chemical  reaction  behavior 
under  non-ideal  conditions 


► 


► 


► 


reactive  potential  not  needed 


multiple  reactions,  multiple  phases  possibl 


reaction  equilibria  not  kinetics 


porous  materials 


Reaction  Ensemble  Monte  Carlo 

Input  requirements: 

•  specify  reactions  occurring 

•  intermolecular  potential  models 

•  ideal  gas  intramolecular  contributions 

(vibrational,  rotational,  electronic) 

predicts  physical  effects  not  chemical  effects 
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cf  cf  cf  cf 
cf  cf  cf  cf 
cf  cf  cf  cf 
&  &  &  & 
cf  cf  cf  cf 
cf  cf  cf  cf 
cf  c#  cf  cf 


initial  state 


shock  Hugoniot  properties 


shock 


8  f  f 


shocked  state 
(E,  P,V) 


(E0,  P0,V0) 


P 


%  shocked 


Hugoniot: 

E-E0='/2(P+P0)(V0-V) 


initial 


V=l/p 


N20  2N 


Zubarev  and  Telegin  (1962) 
Nellis  et  al  (1991) 
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■  RxMC  (non-reactive) 
•  RxMC  (reactive) 
a  experiment 
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a  experiment 
RxMC 
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Schott,  Shaw  and  Johnson  (1985) 


pressure  [GPa] 


Nitromethane 


►  prototypical  EM 

►  27  reacting  species 

►  42  reactions 

'(0 

Q. 

CD 
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C.Mader  (2002) 
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Nitromethane 
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Capabilities  of  RxMC 


►  Assess  accuracy  of  theoretical  models 


►  Predictive  tool  for  energetic 
material  development 


►  Supercritical  fluid  phase 
separation 


energy  storage 


Capabilities  of  RxMC 


►  Predictive  tool  for  other 
non-ideal  environments 


carbon  nanotubes 


Capabilities  of  RxMC 


►  Sol-gel  processing  of 
nanostructured  materials 
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